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Abstract 

In this paper the phonon self energy produced by anharmonicity is calculated using second 
order many body perturbation theory for all bcc, foe and hep transition metals. The symmetry 
properties of the phonon interactions are used to obtain an expression for the self energy as a sum 
over irreducible triplets, very similar to integration in the irreducible part of the Brillouin zone for 
one particle properties. The results obtained for transition metals shows that the lifetime is on 
the order of 10^^'^s. Moreover the Peierls approximation for the imaginary part of the self energy 
is shown to be reasonable for bcc and fee metals. For hep metals we show that the Raman active 
mode decays into a pair of acoustic phonons, their wave vector being located on a surface defined 
by conservation laws. 
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Harmonic phonon calculations based on density functional theory are nowadays rou- 
tinely performed for bulk solids. The dynamical matrix is either obtained from density 
functional perturbation theory^ or from supercell calculationa^i^. To go beyond the har- 
monic approximation quasiharmonic calculations are usually performecl^i^. However in this 
effective theory the phonons do not have lifetime. Ab-initio anharmonic calculations taking 
into account phonon-phonon interactions explicitly are rather rare. There are noticeable 
exceptions with for example the calculations in the diamond structure of Si and Ge^*^, and 
the recent study of graphite by Bonini et al~. Such calculations gives relevant informations 
about the phonon-phonon interactions which may be hidden by the electron-phonon interac- 
tion in experiments. It is important, for example, in the understandig of energy transport in 
thermoelectricity. Looking at the self energy of simple basic elements is therefore of interest. 

In this paper we study the bcc, fee and hep transition metals for the first time. The 
phonon-phonon self energy is calculated for all metals in the crystallographic structure stable 
under normal conditions. The necessary information is then extracted to obtain the decay 
path for a selected phonon in the Brillouin zone. 

The paper is organized as follow. To explain the methodology of our calculation we first 
define irreducible triplets of wave vectors from the symmetry of the phonon-phonon coupling 
function. A formula is then obtained to calculate the (g, u) resolved self energy in a simple 
way. The results of these calculations for the transition metals are then analyzed using 
band decomposition and conservation surfaces for the phonons with the shorter lifetime. 
An approximation proposed by Peierls is also discussed. Finally the phonon decay path 
generating the Raman damping is described for hep metals. 

The strength of the interaction, J-", between phonons of wave vectors q, q', q" in bands 
p,p',p", is given in terms of the eigenvalues, ujp{q), and eigenvectors, e™(g) of the harmonic 
hamiltonian, as well as the third derivative of the potential energy, ^QriR^n Rsts"' 
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In the above equations, are any lattice vectors of a crystal containing cells, and Ar^^ 
a displacement of atom r with mass in cell R in the direction a around the equilibrium 
position. 

According to the second order many body perturbation theory, the third or- 
der hamiltonian Hj, produces the self energy Sp(g, cu) = Ap(g, cj) + iVpi^q^uS) with 



rp(g,w) = J2p'p"'^PP'p"il^^) and 
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The function / is the temperature dependent part given in term of the Bose-Einstein occu- 
pation factor Ugp by 

{pp'p"} 

f = {Uqipi — nqiipii)5{uj + Uqipi — Ugiip") + 1/2(1 + Ugipi + nqiipii)6{uj — UJqipl — Ugiip"). 

{qq'q"} 

However in this paper we are only concerned with the T = hmit where the first term 
vanishes. 

The calculation of the self energy can be greatly improved if we use the symmetry prop- 
erties of the coupling function J^. Let us denote by V the set of permutation operations 
V = {1, -P23, -P12, -Pi3, -P12-P23, -P13-P32} where Pij switch the i and j element of any triplet. 
For example Pi2{a, b, c} = {b, a, c}. The set of rotations of the point group of the crystal is 
called TZ. To make the equations more compact when such rotations are applied to a triplet 
of wave vectors {qq'q"}, we will use the notation R{qq'q"} = {RqRq'Rq"} V-R G TZ. The 
first vector of such a triplet is written as R{q\q'q"}- 

From the definition ([1]) of the coupling function J-", it is then straightforward to show 
the following properties, 

P{pp'p"} {pp'p"} , , 
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Using the invariance of the potential energy under the space group operations of the crystal 

and the law of transformation for the eigenvectorslfi, one can also show that, VP G 7^ 

{pp'p"} {pp'p"} ^ // ^ ..N 

-A = It q + q+ q=G, (5) 
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where G is a reciprocal lattice vector. 

The symmetry properties (E]) and can now be used to define a set of irreducible 
triplets of wave vectors {qq'q"}. By definition a set of irreducible triplet is a minimal set 



of triplets {kk'k"}, which sum up to a reciprocal lattice vector and can be used to generate 
any triplet {qq'q"}, which also sums up to a reciprocal lattice vector, by application of the 
elements of V x TZ. In short {qq'q"} = RP{kk'k"}. 

According to this definition it is sufficient to calculate the coupling function J-" for a set 
of irreducible triplets since all other can be deduced from it. If {qq'q"} = RP{kk'k"}, then 
one has 

{pp'p"} ^P'Hpp'p"} 

T = X . (6) 

{gg'g"} {kk!k!'} 

By analogy to the reduction of Brillouin zone integration to its irreducible part for 

properties such as the density of states, formula ([2]) for the self energy is now reduced to a 

sum over irreducible triplets. Since g" is determined from crystal momentum conservation, 

the summation over g' can be seen as a sum over all triplets starting with wave vector q. 

Then one can write the contribution to Y^ifi^bS) for a phonon decaying to bands p' and p" 

as 
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The ^1 function is different from zero if g = PR{k\k'k"} and is given by the reciprocal of the 
number of times the triplet PR{kk'k"} has been generated by application of the operations 
of P X 7?. on {kk'k"}. The second line is obtained using the symmetry properties Q, (j4j), ([5]) 
and RP = PR. The weight coeficients Cp are calculated once for all from Cp{q, {kk'k"}) = 
^^g^ 5i(g, P-R{A;|A;'A;"}), and it is also usefull to remark that Ci = Cp^^, Cp^^ = Cp^.^p^^ and 
Cpj3 = Cp^2P23- Equation ([HD is particularly useful for computer calculations since it can 
easily be parallelized over irreducible triplets of wave vectors, as with single /^-point for one 
particle properties. A computer code has been implemented from these equations, and in 
the following it is applied to bcc, fee and hep transition metals. 

With the exception of manganese, the transition metals crystalize in the bcc, fee and hep 
structures. The second- and third-order force constants can be seen as derivatives of potential 



energy or dervatives of forces which are obtained from these structures using first principle 
calculations. In particular the third order force constants are third derivative of potential 
energy with respect to atomic displacements. Therefore they are calculated from forces on 
atoms in a supercell containing two atomic displacements. The total number of atomic 
displacement pairs is reduced using crystal symmetry. In our study, we employed finite 
displacement method to calculate the derivatives, but to improve the accuracy displacements 
of plus and minus directions are applied if they are not symmetrically equivalent. The third 
order force constants are usually overdetermined in this way. The tensor elements are then 
determined using pseudo inverse, which is the technique also employed for second-order force 
constants^. The details of these calculations are given in the appendix. 

To obtain the electronic structure and forces we employed the projector augmented wave 
method^i, in the framework of density functional theory, within the generalized gradient 
approximation of Perdew, Burke and Ernzerhof^ as implemented in the VASP code^^"— . 
Spin polarized calculations are performed for Fe, Cr, Ni and Co. The supercells of the 
bcc, fee and hep structures contain 16, 32, 16 atoms respectively, and are limited by our 
computational resources. A plane- wave energy cutoff of 300 eV is used, and k-point sampling 
meshes of 12 x 12 x 12, 12 x 12 x 12 and 16 x 16 x 8 are used for bcc, fee and hep supercells, 
respectively. The Methfessel-Paxton scheme^fi is employed with a smearing of 0.2 eV. The 
cell parameters are relaxed using until the stresses becomes less than 10~^ GPa. Atomic 
forces are obtained with an energy convergence criterion of 10~^ eV. For some metals, such as 
Cr, the electronic ground state we obtain from density functional theory can be questionable. 
However the forces extracted from those calculations may still be used to calculate the forces 
constants. For example for Cr we have checked that the second order forces constants gives 
a phonons spectrum, at the point N of the Brillouin zone, at most different by 6 % of the 
experimental values^. Since the third order force constants is even more short ranged, we 
assume this approximation to be still acceptable. 

According to equation ([8]) and the force constants previously calculated we can obtain 
the damping functions. They are calculated for all bands and all high symmetry points^^ 
in the first Brillouin zone. These functions are non zero between and 2umax but we found 
the center of gravity located at about 2u}, where u is the average phonon frequency over 
the Brillouin zone. Those are quite smooth functions for bcc and fee metals whereas they 
exhibit a more complicated structure for hep metals. 
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FIG. 1: Imaginary part of the self energy in THz. The values are given for all bands and all high 
symetry points of the Brillouin zones. 



The probability decays of harmonic phonons are then found according to the equation 
l/2Tqp = ImLqp{ujqp). They are represented on figure [1] for all bands and all high symmetry 
points of the Brillouin zones. One can see that the minimum lifetime tend to decrease 
toward the right end of the series. But even if most of the calculated values are greater than 
0.5 X 10"^'^ s, one cannot distinguish a clear trend in their distribution for a given crystal 
symmetry. 

The imaginary parts of the self energy corresponding to the minimum lifetime are rep- 
resented in figure [2] as functions of frequency. The vertical line shows the frequency of the 
irreducible representation it belongs to. They are always located within the lower tail of the 
self energy. 

To better understand how a harmonic phonon of frequency ojq^p^, which will be the phonons 
with minimum lifetime, acquires a finite lifetime qQpo{<^qopQ)i one should remember that this 
quantity is constructed from two parts. The one with the delta functions gives the decay 

{PP'P"] 



processes which are allowed by the conservation laws and the 



gives the proba- 



{qq'q"} 

bility for such decays to happen. The two conservation laws, for energy Uq^^p^ = Uqfp' + Uq"p", 
and momentum qo = q' + q" + G, are coupled equations which define a conservation surface 
in reciprocal space: a phonon in mode goPo will decay in two phonons of bands p' and p" 
with wave vectors having their extremities on that surface. For each metal at least one 
couple of bands has a large probabihty decay. The conservation surfaces corresponding to 
the strongest ones are plotted in figure |3] and the percentages for such decays are given in 
the caption. In such a way one obtains a very clear view of the processes which generates the 
lifetime since we know the bands to which the phonons decay as well as their wave vectors 
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FIG. 2: The imaginary part of the self energy calculated at the point q and band p is represented as 
a functions of frequency. The continuous line is the full calculation and the dashed one correspond 
to the approximation due to Peierls. For the hep metals the functions for the Raman active modes 
are shown with dash-dotted line and their approximation as dotted line. 

(the band indices are given in figure [2]) . For given p' and p" we should however remember 
that when p' ^ p", the surface is always composed of two sheets. One centered at the origin 
where q' is located, and the same shifted by go where q" lies in. For clarity only the first one 
is represented on figure [31 All surfaces we found are open surfaces, and with the exception 
of Ti, they have a tube shape along go- 
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FIG. 3: Conservation surfaces for the bands with the stronger probabihty decays, qo is chosen as 
the vertical direction. The surfaces are generated using the Xcrysden software^^. The band indices 
and the contribution to the damping function (in percentages) at the harmonic frequency are as 



follow : Sc{p',p" = 4,1; 36%), Sc-Raman(p',y' = 2,1; 76%) , Ti {p',p" = 3,1; 27%), Ti-Raman 
{p',p" = 2,1; 50%), V {p',p" = 2,1; 80%), Cr {p',p" = 1,1; 50%), Fe {p',p" = 2,1; 41%), Co 
{p',p" = 4,2; 28%) Co-Raman {p',p" = 2,2; 47%), Ni {p',p" = 2,1; 59%), Cu {p',p" = 2,1; 51%). 



The three surfaces at the second row correspond to Raman decay. 

Now if we consider the frequency ojqopo as a variable parameter u, we generate a family 
of surfaces S{oj) whose shape and area give the joint density of states. 



As an approximation, Peierls^ proposed the damping function Tg^p^^ (u) to be proportional 
to the joint density of states. In fact one can also simply fix the proportionality constant 



where we defined A as Er^r^r-, Ea,a,a^ ER,R,i^or"R2r2,R3r-f ■ I* represents an aver- 



aged measure of the anharmonicity. Uq is the average value of frequencies at point q. 
The approximation is clearly good for fee and even bcc transition metals where the anhar- 
monicity seems to be describable by a single number a for these processes. However the 
approximation is not accurate for hep metals. The coupling function has a much stronger 
dependance on the phonon modes. This comes from the hexagonal structure which has two 
atoms per cell. This combination gives third order force constants which are anisotropic, 
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by its average value and rescaling by oo, 
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and complicated interferences between phonons eigenvectors, accoustical and optical, which 
cannot be removed from the calculation. We have also calculated the Raman damping func- 
tions of the hep transition metals and the same conclusion is obtained. For this structure 
it is the mode which is Raman active. These phonons have a much longer lifetime with 
2.1 X 10~^s for Sc, 5.8 x 10~^s for Ti and 1.8 x 10~^s for Co. Their conservation surfaces 
are presented at the second row of figure [3] and the damping functions are shown in figure 
[2J For this mode the stronger decay to a couple of bands p'p" is much more selective than 
in the case of the modes with minimum lifetime. Remarkably, an optical phonon in the 
Raman active mode will decay into a pair of acoustic phonons in almost all cases. The 
wave vectors of such phonons are located on the surfaces shown in figure [31 These surfaces 
are closed, with very simple shapes. This seems to indicate that simpler models could be 
constructed for these decay processes. 

In conclusion, we have calculated the phonon-phonon self energy of bcc, fee and hep 
transition metals. The decays for the phonons with minimum lifetime were studied and the 
conservation surfaces calculated. We found that for bcc and fee metals the imaginary part 
of the self energy is approximately proportional to the joint density of states whereas this 
approximation fails in the case of hep metals. The Raman damping was also examined for 
these metals and we found that a phonon decay into a pair of acoustic phonons whose wave 
vectors are located on spherical-like surfaces. 

The authors gratefully acknowledge the French Agence Nationale de la Recherche (ANR) 
for financial support under contract # 07-MAPR-0015-04 as well as the Grants-in-Aid for 
Scientific Research (A), Scientific Re- search on Priority Areas (Grant No. 474), and the 
Global COE Program, all from MEXT, Japan, and the PIE - Programme Interdisciplinaire 
Energie of CNRS, France. 

Appendix A: Computations of second- and third-order force constants 

The potential energy of a phonon system is represented as a function of atomic positions, 
^(r_Riri5 • • • 5 ^Rnt„)j where r^j,- is the atomic position, and n and are the number of atoms 
in a unit cell and the number of unit cells, respectively, and Ri are the indices of atoms 
in a unit cell and the indices of unit cells. 

A force on an atom is the first derivative of the potential energy with respect to an atomic 
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position, 

dV 

Fnr = (Al) 

a, P, are used for the indices of Cartesian coordinates. A second-order force constant 
<l>"'^ is the second derivative of the potential energy as function of atomic positions, 

and a third-order force constant ^"^"^ is the third derivative of the potential energy as 
function of atomic positions. 
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Using finite differences, the derivatives in Eqs. (lA2p and (]A3p are approximated by 



^RiTi,R2T2 - .8 ' y-^^) 



R2T2 

and 



a/37 ^ ^^flin.flarat'^^jiaraJ 



RsT3 

_ ^RiTi,R2T2i^' Rsrai ^RiTi,R2T2 

respectively. Ar^ and Ar'^ correspond to the finite atomic displacements. The Ar'^ and Ar'^ 
appering in the parentheses of forces and force constants mean that the values are calculated 
under the displacements. 

To compute the second-order force constants, we employed the technique presented by 
Parlinski et. al^ and the third-order force constants are obtained in a similar manner. In 
the following sections, the computational details are given. 
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1. Computation of second-order force constants 



Second-order force constants are computed through the approximation (IA4p with small 
displacements. For computational convenience, a second-order force constant tensor for a 
pair of atoms, RiTi and R2T2, and an atomic displacement are represented by a 9 x 1 matrix 
P and a 3 X 9 matrix U given by 



P(i?2'r2, -RiTi) = [$^^ $2^^ ^yy $2^^ <|)^^ ^^y $^ 



(A5) 



and 



'^l 0^ 



U(i?2r2) 



(g) [/^^x ^^y 



(A6) 
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010 

respectively. Using these matrices, a force on an atom, which is in the form of a 1 x 3 matrix 
F, is obtained by 

F{R,n) = -U(i?2r2)P(i?iri, i?2r2). (A7) 

Simultaneous equations of different atomic displacements for a pair of atoms are then com- 
bined as 



V ■ J 



U2 



(A8) 



With sufficient number of atomic displacements, Eq. f lASp may be solved by pseudo inverse 
such as 

U2 



/ttA Vf,\ 



V '■ J 



(A9) 



However with the help of site-point symmetry, the required number of atomic displacements 
to solve the simultaneous equations may be reduced. If -R'it( is the image of atom RiTi by 
a site-point symmetry operation of atom R2T2, Eq. (IA7p becomes 



F(i?;r- 



-\J{R2T2)V{R[t[,R2T2) 
-\]{R2T2)AP{R^n,R2T2) 



(AlO) 
(All) 



where F(i?'^r() is the force at the atomic site obtained from the original atomic site by the 
site-point symmetry operation, and A is the 9x9 matrix that is used to rotate P along the 
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site-point symmetry operation. Using Eq. (lAlip . the combined simultaneous equations are 
built such as 



(2) 



^(2) 



UiA(2) 

U2A(2) 



v 



(A12) 



V ■■ } 

where the superscript with parenthesis gives the symmetry operation index. This is solved 
like Eq. flX9|) . 



2. Computatation of third-order force constants 

The finite difference approximation for the third-order force constants is represented by 
matrices as 

AP(i?iri, R^T2) = V{R^T3) ■ Q{Rin, R2T2, R3T3) (A13) 

where AP, V, and Q are the 9 x 1, 9 x 27 and 27 x 1 matrices corresponding to A^"'', Ar'^, 
and respectively, and are given by 



AP^+3(a-i) = A*-/^ a, (3 = 1,2,3 



(A14) 



and 



Q^+3(/3-l)+9(a-l) = ^'"^^ a,(3,-f= 1, 2, 3 

fl o\ 



(A15) 
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(A16) 
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respectively. Simultaneous equations are constructed by a similar manner to Eq. ( ]A8|) as 



APi 
AP2 

V • / 



V2 
V • / 



Q 



(A17) 
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This may solved by pseudo inverse such as 



Q 



V2 



AP2 



(A18) 



The number of pair of displacements to calculate can be reduced using symmetry operations 
that conserve a third-order force constant tensor for a triplet of atoms. If -R'it{ is the image 
of atom RiTi through a symmetry of the displaced structure, then one has 



AP(i?;r(,i?2r2) 



V(i?3r3)-Q(i?'in,i?2r2,i?3r3) 
V(i?3T3)-B-Q(i?iri,i?2r2,i?3r3) 



where B is the 27 x 27 symmetry operation matrix that transform the tensor Q. 
The simultaneous equations are then written in a similar manner to Eq. (]A12p as 



APi' 



(2) 



AP 
AP 



(1) 
2 

(2) 



ViB(2) 
V2B(i) 

V2B(2) 



Q 



(A19) 



and this is solved in the same way as Eq. (lAlSp using the pseudo inverse method 
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